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, Thin viscous Keplerian accretion disks are considered asymptotically stable, even though they 

£\j . can show significant dynamic activity on short timescales. In this paper the dynamics of non- 

axisymmetric hydrodynamical disturbances of disks are investigated analytically building upon the 
steady state three-dimensional structure and evolution of axisymmetric perturbations explored in 
previous work. Assuming a polytropic equation of state solutions are found by means of an asymp- 
totic expansion in the small parameter measuring the ratio of the disk thickness to characteristic 
radius. In-depth analysis shows that every perturbation that disturbs the radial velocity induces 
significant transient growth in the (acoustic) energy of the evolving disturbance. This effect is most 
evident in the density and vertical velocity. The transient growth observed is tied to the non- 
[•j - ] ' separable nature of the solutions where, in particular, pattern evolution is controlled by a similarity 

, variable composed of the radial coordinate and time. This leads to growing winding perturbations 

HH ■ that display successive radial peaks and troughs. We argue that these transient non-axisymmetric 

pH ' structures may precipitate secondary instabilities which, consequently, may be a critical element for 

Q a new alternative picture of turbulence arousal in non-magnetized astrophysical disks. 

O ' PACS numbers: Valid PACS appear here 

& 

CO 

! I. INTRODUCTION 

Since the late 1980s a new perspective has developed in hydro-dynamic stability theory. This new paradigm arose 
from the long-standing problem of linearly stable shear flows that experimentally exhibit transition into turbulence. 
The conventional approach had been to examine the linear stability of fluid systems via modal analysis ("normal- 
modes"). In practice this means that the determination of the eigenvalues and eigenfunctions of the linearized 
perturbation equations of a given flow indicates the time asymptotic behavior of the disturbances and, consequently, 
helps to determine the long-time stability of the base flow. One of the best expositions of this classical approach is 
found in the book by Drazin & Reid Q. 

The origin of the new perspective can be traced to the fact that linear stability analysis of shear flows in general gives 
rise to non-normal operators, i.e., linear operators that do not commute with their adjoint. A typical modal analysis 
of a problem governed by a non-normal operator can lead to an incomplete description of the full breadth of responses 
■ possible for the original initial-value problem (IVP). For example, non-normal operators will have eigenfunctions that 
are non-orthogonal, and/or imply the existence of solutions which are unobtainable analytically (see the discussion 
in Chapter 8). 

Especially critical in this matter is the fact that the non-orthogonality of the eigenmodes can lead to transient 
dynamics which, in turn, can imply strong transient growth (TG) for suitable initial conditions, e.g., in perturbation 
energy or enstrophy. The implications of this for a variety of shear flows has been studied in numerous earlier works, 
e -g-j S EL HI- The framework for stability calculations has thus shifted from just focusing on the time-asymptotic 
behavior of a disturbance towards studying its TG as well. The recent review article by Schmid [|| gives an up-to-date 
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account on non-modal stability theory and its successes. A detailed exposition of the subject, including its different 
aspects and possible extensions, can be found in the book by Schmid & Hennigson Q- 

Ioannaou & Kakouris Q were first to apply this non-modal perspective for an astrophysical setting, namely, for the 
problem of accretion disks ( "ADs" and "AD" for singular usage) . Accretion disks are both important and ubiquitous 
astrophysical objects which are thought to power systems as diverse as young stellar objects, close binary systems 
and active galactic nuclei. ADs are flattened, swirling flows in a gravitational field of a central compact object, 
wherein high specific angular momentum fluid accretes onto the central object. In order to reconcile theoretical 
models with observations, an efficient mechanism is needed that dissipates energy and transports angular momentum 
because the slow spiraling-in and eventual accretion of fluid depends critically on this process. When "viscous" ADs 
were theoretically proposed [ij, [To| it was recognized that an anomalous dissipation and transport mechanism must be 
present in ADs since their hydrodynamical Reynolds numbers (Re) are enormous. Because fluid turbulence greatly 
enhances transport, turbulence has been proposed to operate in ADs as angular momentum can be transported in 
rotating flows with the help of a turbulent eddy viscosity. To date, a detailed theoretical understanding of turbulence 
and the transition to it is still lacking. Consequently, the effective viscosity in disks has been approached in a 
phcnomcnological way through parameterizing the effective viscosity coefficient with the help of a non-dimensional 
paramcter(o;) on the basis of dimensional arguments. This simple approach has been exceptionally fruitful, giving 
rise to successful interpretations of many basic observational results (ill . [l2| . 

Robust hydrodynamical stability criteria like the Rayleigh and Solberg-Hoiland criteria indicate that thin non- 
magnetized Keplerian ADs are linearly stable. To date there are no demonstrations of long-time dynamical activity 
in global hydrodynamical simulations of AD flows with sufficiently high Re. On the other hand, high resolution 
3D numerical simulations of local disk sections [l3[ report that a subcritical transition to turbulence does exist at 
very high Reynolds numbers despite suggestions that Coriolis effects quench such a transitionflil fl5|. However, the 
efficiency of turbulent transport in such subcritical flows appears to be insufficient to explain the transport implied 
by the observations of ADs. In this study we use the term global when the calculation includes a sizable portion of 
the AD. In contrast, a local calculation is one that is based on the so-called shearing box approximation fTolll7l.[l8| in 
which calculations are performed on small "Cartesionized" sections of the disk. 

In the years following the work of Ioannaou & Kakouris 0] (who had employed a global approach albeit to a two- 
dimensional configuration) a number of research groups have used the local approximation to study the relevant IVP 
in various settings [TH, 0, H| HH, H2, HI, HH Hf| . A common conclusion reached by these studies is that TG may be 
copious and that, under the right conditions, nonlinear interactions may give rise to what is called a bypass transition 
to turbulence. It is worth mentioning that vortices and spiral waves are singled out in some of these works as being 
instrumental for such a transition. 

This work is intended to further understanding what the non-modal IVP perspective can teach us about the global 
transient dynamics of thin slightly viscous hydrodynamical Keplerian disks. To be specific, we consider flows with a 
viscosity (perhaps of turbulent origin [l3|]) that is insufficient to provide the angular momentum transport implied by 
observations. We are interested in examining the excitation of global secondary flows that can, in turn, possibly give 
rise to angular momentum transport and/or create conditions for a secondary instability atop the weakly turbulent 
state [II, [23, [II, [II [II. 

In the context of ADs one is confronted with a vast system on which there is no experimental control. In lieu of 
this, Ionnaou & Kakouris Q reported in their study that stochastic forcing was found to lead to persistent activity 
with angular momentum transported outward. More recently, Zhuravlev & Shakura [3l| applied optimal perturbation 
strategy to two-dimensional sub-Keplerian toroidal configurations. They found that optimal perturbations giving rise 
to substantial TG are composed of certain combinations of non-axisymmetric eigenmodes. 

Our tactic is to study the transient dynamics of specific three-dimensional perturbations of a Keplerian AD in which 
the vertical structure is taken into account. The base flow is the Kluzniak-Kita [lH (hereafter KK) analytical solution 
of a polytropic steady viscous axisymmetric disk. This solution was obtained by representing all dependent functions 
as an asymptotic series in the small parameter e which measures the disk's "thinness" (i.e., the ratio of the disk's 
vertical thickness to its typical radial scale). The lowest order solution is the classical Shakura-Sunyaev solution [j|, 
while higher order terms provide the velocities and structure functions in the meridional plane, as well as corrections 
to the Keplerian angular velocity [33j. This kind of perturbation strategy for ADs had first been introduced in the 
context of AD inner boundary layer [34j. We wish to reiterate here that implicit in the KK solution is that there is 
some viscosity on the smallest scales, perhaps due to the subcritical transition discussed previously [HI]. 

Umurhan ct al. [35[ (hereafter UNRS), extended the analytical KK solution to consider its time dependent ax- 
isymmetric response. The analysis of the corresponding IVP indicated that typical initial data naturally exhibits 
transient growth lasting many rotation periods before ultimately decaying. The lifespan of this response is inversely 
proportional to the viscosity parameter a, which is consistent for studies of simpler systems [3]. Because UNRS 
examined only axisymmetric initial disturbances and flows, the transiently growing patterns could not induce any 
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effective radial angular momentum transfer. In addition, no obvious mechanism was identified that could serve as a 
suitable secondary instability candidate. 

In this paper we report on a work based on the same ideas applied to the same slightly viscous system, but allowing 
for non-axisymmetric disturbances. Consequently the results, which remain analytical as in UNRS, exhibit a far richer 
structure. The asymptotic expansion procedure used in UNRS is applied in which the same small parameter (i.e., the 
disk's 'thinness') is exploited resulting in an IVP that is analytically treatable. The expansion procedure employed 
here develops a finite-amplitude nonlinear solution. Thus, although the higher order terms and time-dependencies arc 
governed by linear operators, the solutions developed are finite-amplitude and they should not be confused as being 
infinitesimal solutions. 

The use of approximation methods and the simplistic polytropic assumption have the obvious advantages that 
the treatment can be analytical and the responsible physical effects leading to any interesting dynamics may be 
transparently traced. It is clear, however, that the present analytic analysis ultimately should be complemented 
with a detailed and uncompromising numerical solution which includes a proper treatment of energy generation and 
transfer (36J. We also stress that we do not aim to obtain a general solution for the IVP, nor do we intend to itemize 
all possible solutions. We also do not aim to find the optimal perturbation, which would be a strategy appropriate 
for a laboratory transition study. The purpose of this study is to show that there exists a solution that gives rise 
to prominent transient growth in the energy when small but non-infinitesimal disturbances are introduced. We also 
investigate how the different parameters of the system (i.e., a, symmetry) affect this growth. Implicit in our approach 
is the tacit assumption that the possible perturbation spectrum in a realistic AD is so rich as to allow essentially any 
initial condition we desire. We are specifically interested in finding a global transiently growing spatial pattern of the 
density and assessing its physical consequences. 

The paper is organized in the following way. In Section II we formulate the problem by stating the assumptions, 
notation and scaling, leading to the basic set of non-dimensional partial differential equations for the flow. In Section 
III the asymptotic expansions of the various dependent variables are given and the equations are solved in each order 
in the small parameter e, up to the second order. As already stated, the lowest order gives the steady Shakura-Sunyaev 
solution and the higher orders provide analytically the time dependent solution for a given initial disturbance. The 
most interesting effect — the transient growth in the density and vertical velocity — occurs at the second order (Section 
III-D). In Section IV and V, the solutions are discussed in some detail and we examine with the help of graphic 
visualization some of their important physical properties and finally, in Section VI, we summarize the results and 
their meaning. Since the analytical procedure needed for obtaining the solutions, in various orders, is lengthy and 
involved we defer some technical details to the Appendices. 

II. FORMULATION OF THE PROBLEM 
A. Notation, scaling and basic assumptions. 

The problem is formulated in cylindrical coordinates (r, z, <f) with P and p denoting the pressure and density 
functions. The cylindrical components of the velocity are u, v and rfi, where f2 is the angular velocity. Additionally, 
c s is the sound speed and rj the dynamic viscosity. The accretion flow is in the Newtonian gravitational field of a 
central object with point mass M. 

The coordinate r is scaled by its value at a typical point r* (we shall denote the scaling dimensional variables 
by an asterisk) and the density by a typical density value p*. The first assumption is a polytropic pressure-density 
relation P = Kp 1+1 ^ n , with n being the polytropic index and K a constant. Consequently the pressure is scaled 
by P* = Kpl +1/n . This also gives the typical sound speed (squared) c 2 sif = dP*/dp*. We shall scale the meridional 
velocities with c st and express the rotational angular velocity in units of its Keplerian value at r*, so that fl 2 = 
GMjr\. This then allows us to express the vertical coordinate scale, h*, using the second assumption, that the flow is 
approximately in vertical equilibrium due to thermal pressure support (that is, ft* = c s */f2*). Our third assumption 
is that the disk flow is cold, or, equivalcntly, that the azimuthal rotational velocity is highly supersonic. This is 
equivalent to assuming very efficient cooling, by , e.g., radiative losses from the surface of the disk. As a result, we 
find that e = ft*/r» = c s */(il*r») <C 1, which happens to also measure the disk's thinness, as we alluded to in the 
Introduction. This small number e will thus appear in the non-dimensional equations characterizing the flow and 
becomes our natural expansion parameter. 
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B. Equations 



The non-dimensional polytropic relation is P = p l+l / n and we define, for convenience, a function W = j dP/ p. 
This gives supplementary polytropic relations, which will ultimately allow to express P, c s and W in terms of the 
density alone. 

( 1 ) 



: n^- = nc 2 = (n + l)p V ™ because c 2 s = (1 + -)p 1/r 
dp n 



The full (time-dependent, non-axisymmctric) scaled hydrodynamic equations in cylindrical co-ordinates read 



' i \ r _ „ , „ , 2 



-1 (l + e 2 ^ 



I - J {d z (r)d z u)} 



I — i <j -pd r W + d z (rjd r v) - ^d r {r)d z v) + dj, (r)d r Q) - y <9«^ - ^d r 



^^ x - 



Jfryu + 2^ ^ r Q rU ^ _ 2^ 1d r {ru) + \d<p (rjd^u) 
p J [ r z r 6 Lr J r z 



a t o + vd z n + no^n = a, (»j5,n) - e {^<9 r (r 2 o)} + 



pr 2 " ) { \ dr i^drty - pd^VF + |efy (t7^S1) + <9 2 (77^ u) - ^ (rjd z v) [• -+ 



1 

pr 2 



12 4 

-<9 r (rrjdfu) - -d<f, (r]d r u) + — fy, (tju) 

7 o o/* 



(2) 



(3) 







2 


1 






1- - 




r 2 ) 







4 2 

-<9 2 (r?<9 2 u) + (?7<9 2 f2) - -<9 2 {■qd^Vt) 
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dtp + d z (pv) + 8$ (pfi.) + e-d r (rpu) = 0, 

r 



(5) 



where (J2H|) are, respectively, the radial, azimuthal and vertical momentum conservation equations, while (O is the 
equation of mass conservation. Note that this set of equations is very similar to the one in KK and UNRS. The 
terms containing angle derivatives are missing in KK and UNRS, because of the axisymmetry assumed there, while 
the terms containing time derivatives are missing in KK, who considered only a steady state. 

We impose a Lagrangian pressure condition on the vertical surfaces (namely, that the pressure on the surface be 
zero). We also require that the stresses on the vertical surfaces be zero as well. See Appendices C-D for details. The 
radial boundary conditions require some discussion. We will consider these equations for an extended ring, in which 
the internal radius is considerably larger than the zero-torque radius of the disk r+, that is, r* ^> r + (sec KK for a more 
extended discussion) and the external radius is significantly smaller than the disk outer edge. In this way we can avoid 
the treatment of inner and outer boundaries, which greatly complicates the problem, presumably without changing 
the substantial result for the bulk of the disk. Still it would be interesting to address in the future the dynamics at the 
inner edge region in the case of black holes (using a different scaling and possibly matched asymptotic expansion) and 
in the boundary layer in the case of neutron stars (by means of matched asymptotic expansion). The outer boundary, 
at which the disk is fed by mass, depends strongly on the astrophysical system in question and is of interest as well, 
but we defer also this question to future work. 



III. ASYMPTOTIC SOLUTIONS 



A. General 



The asymptotic expansion approach consists of expanding all functions in powers of e (e.g., KK, UNRS) as follows 

f(r,z,<j>,t) = fo(r,z) + e fx(r,z,<j>,t) + e 2 f 2 (r, z,<j>,t)... (6) 
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We postulate, similarly to UNRS, that the expansions of all the dependent variables are such that the terms fj can 
be split into a spatial steady part fj — the steady base flow — and a time dependent dynamical disturbance /j 

fj(r, z, 0, t) = fj(r, z,<j>) + fj(r, z, <f>, t) (7) 

This form presupposes therefore that the time dependence is included only as an additive function and from the first 
order in e and on. We remark here at the outset that the viscosity, rj, is not expanded. Instead, we shall express it, 
using a standard prescription (as found in the Shakura-Sunyaev solution 0]) by the lowest order dependent variables 
(see section III C). 

Without assuming anything additional we obtain what we call the complete set of equations in the first three orders 
in e — see in Appends. IB 11 - IB 31 respectively. We then take the steady base flow to be the KK solution and so it is 
axisymmetric. Although, this is very similar to the UNRS approach, we shall not exclude a priori from the expansions 
terms which were excluded by UNRS, who set them identically to zero, in accord with the work of KK. These include 
the following variables: u 0l Qi, vq, v±, p±, and thus Pi and W\. Along our work we shall explicitly state which 
assumptions of this kind can be "derived" and which arc a matter of choice and thus are assumed. Obviously, this 
work differs from UNRS also by allowing non-axisymmetric disturbances. Note that the vertical equilibrium at the 
zeroth order followed from the axisymmetry of the base flow — explicitly of po- 

Before turning to the equations and their solutions in the three lowest orders in e we assume, for the sake of 
simplicity, that n = 3/2 (as in KK). For our calculations we allow for general values of n and we agree with UNRS 
who find that the results are little influenced for reasonable n values. For n = 3/2, as we assume henceforth, we get 
from the non-dimensional polytropic relation and equation ([1]) that 

D 5/3 T ,r 5 2/3 2 ^ 2/3 , c n 

Po = P ' , W = -p > , c s0 = -p < . (8) 

Finally, we remark that if g(p) is a smooth function of the density alone, and the density is written as p = po + Sp, 
where Sp is a small perturbation atop po, then the perturbation in g, i.e., Sg is well approximated by Sg = (dg/dp)o5p. 
Thus 

This result will be used in the orders e and e 2 below. 



B. Order e 



The complete equations at this order are simple because there is no time dependence at this order. Since the 
steady base flow is also axisymmetric, it further significantly simplifies the equations and gives a particularly compact 
equation set. 

nl = 1 (io) 

dn d ( dn \ 
V0P °-dz- = Tz\*W) (n) 

8v z dW 4 1 d ( dv Q \ 

V 01T- = 3 a ^ q t - ( ) ( 12 ) 

oz V s oz 6 po oz \ oz J 

d(povo) = Q 
dz 



(13) 



Solution 



Equation (jTUJ) guarantees the Keplerian form f2 = r 3 / 2 , which makes equation (fTl) trivial, and equation (TT 
immediately gives 



PoVq = f(r)- 
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Since vo is an odd function in z (and po is even in z), we get that vq = identically and equation (p~2|) can be easily 
solved (subject to the boundary condition po(r, h) = 0, where h(r) is the height of the disk) to yield Wq(t, z) and thus 
Po(r, z) and po(r, z). Thus the O (0) solution is identical to the KK solution at the same order, 

^0 = r-*l\ v Q = 0, p (r, z) = (^f) ' , (14) 



with h = h(r) as given in KK and in UNRS, 

1/6 



kij) = A (l- A /I± ) with A 



r \ V T 



M / 80 /5 
~o~ I 37rV 3 



-,1/6 



(15) 



Also 



following from 



^ 5/2 b(h 2 -z 2 



C. Order e 1 

The complete (e) system, as given in Append. IB 21 is significantly simplified, when the axisymmctry of Sl an d 
the fact that «o = 0, from the zero-order solution pT|) . are used. 

Taking only the " unperturbed" parts gives rise to the following time independent equation set 

LI. ( fb*\ 



-2r£l Q Sli = ——[r,^\ (17) 
po dz \ dz J 

r 2 dr po dz \ dz J 

am + j_d_ ( dvA 2^ ( + J_9_ (duo) (19) 

dz 3po dz \ dz J 3 por dz \ dr J p r dr \ dz J 
ld(rp u ) d(povi) 



dr dz 



= (20) 



These equations are for the axisymmetric base flow (identical to KK; UNRS did not consider this order and took all 
the contributions to be zero). In any case, subtracting this from the complete equation set gives 

- 2rfi fii = (21) 

This equation guarantees that = and we can use the notation fii for f2i, because it is time independent in 
this case. The resulting second equation gives no additional information (time-dependent, and see below) and the 
remaining time-dependent equations to this order are 

dv[ dv[ dw{ , a a ( dv[\ 

= —aV + 3p- dz ) (22) 

dp'i , o <Vi djpovl) 



We note the appearance of the operator 



dt d(j) dz 



D * S ! + n 4' (24) 



which we shall use henceforth, for economy of notation, remembering also that f2o( r ) = r 3//2 - The subscript "0t" in 
the definition of this operator has the purpose of reminding that it contains both the time and angular derivatives. 
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1. Solution 



In the previous order the explicit form of the viscosity rj was not needed, but now the situation is not as comfortable. 
Following Shakura and SunyaevQ, we posit the following form of the (cf>, r) viscous stress tensor component 

M=aP , (25) 

where a is an adimcnsional parameter. This gives 

= aP Q -» r) = ^ar 3 /' 2 p 5 /3 (26) 

Notice that this prescription is vertically dependent, i.e., the vertical distribution of stress is assumed to follow 
the pressure distribution. We now substitute this r/ together with the polytropic relation, f^o = r~ 3 / 2 and use 
d z Wo = —z/r 3 (from equation ([12"]) ) in equations (|17H18[) . which become, after some straightforward algebra 



dfl 



dr 



^/ V ^_!^ + (2ri5l) = (27, 

3 oz z 6 Oz 

Now, since we consider only even (in z) solutions for uq and fii, it can be shown (see Appendix [Aj that the only 
bounded solutions are the trivial solutions: 

uo(r,z) = Qi(r,z) = 0. (29) 
Equations (|17til8[) thus become trivial and equations (|19H20[) assume the form 

= -^ + A#(V^ (30) 



dz 3po dz V dz 



= (31) 

oz 

Considering equation (|31[) first, we get poi'i = f(r), but since v\ is an odd function of z (i.e., it is zero at z = 0), this 
implies v\ = identically. Similarly, from equation (|30[) we get W\ = f{r), but at the disk vertical edge (some large 
enough z) this function must be zero. Consequently we will have, similarly to UNRS, 

v 1 (r,z) = Wi(r,z) = Pl (r,z) = 0. (32) 



We are thus left with the equations for the perturbations (|22[) - (|23[) . UNRS omitted in their expansion the first order 
terms of the vertical velocity and density, as well as their time dependent perturbations. Here we shall allow for these 
perturbations and therefore will have to solve equations (|2"2")l - ([2"B")) . These equations can be rewritten as 

t> ' 5 d ( -V3 A , 4 9 / fli/A 

^^ = -3^l Po N + toTzV-te) (33) 



Vfa Pi = t, . (34) 

oz 

]/3 ______ 

where we have used W{ = (2/3)(Wo/ po)p[ = (5/3)/o ' (see equation ©). Applying again V^ t to and 
substituting (|3^|) we get the single equation 

^ / 4 a / a<\ s 2/3 a 2 , 25dUf) a , s a 2 ( P g /3 ) , 

^ " 3^ *^ ^ arj - 3 Po ^ ~ -t-BTa? 1 ~ = °' (35) 

where r\ is given in terms of po by equation (|26p . The above differential equation can be written as 

Cv[ = 0, (36) 
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where the form of the linear differential operator C can be easily inferred from the above. This form becomes explicit 
using equation ([T4"|) 

P0 (r,z) = (5r 3 )- 3 / 2 (/ l 2 -z 2 ) 3 / 2 , (37) 
the viscosity in terms of the coordinates (and constants) 

rj(r, z) = \ a r^ P l" = J -^ r2 ar-Hh* - z 2 f' 2 , (38) 

and remembering that h is a known function of r (see KK, UNRS). 

Equation (f36|) can be solved analytically. At this purpose we introduce the similarity variable T = flat = t/r 3 / 2 
and the Ansatz v[ — £>i(r, z, <fi)e sT + c.c, where s is an eigenvalue. Note that v[ must be periodic in 0, with period 
27r, and thus it can be Fourier expanded and the Ansatz takes the form 

oo 

v[<?,z,4> > t)= E ii( m )(r,Oe sT+lm * + c.c, (39) 



m— — oo 



where £ = z/h(r) is an another similarity variable and the lower index in parentheses denotes the Fourier components. 
Now, for each particular Fourier component, the Ansatz and the resulting eigenvalue equation are like the ones 
discussed in UNRS (see equation (48) of that paper, which is for the case of axisymmetry, i.e., m — 0). The procedure 
is thus identical. 

We summarize here the main result (which will also be used in the discussion of the next order). The ordinary 
differential equation resulting from the substitution of the above Ansatz, is known as the Gcgenbauer (or hyper- 
spherical) equation and its solutions are known in terms of combinations of the associated Legendre functions (also 
known as Gegenbauer polynomials), e.g., H3, [H|]. These analytical expressions are also included in Wolfram's 
Mathematica 6 software, which makes our calculations feasible. Applying the upper (and lower) boundary conditions 
gives a condition on the eigenvalues. All the eigenvalues thus obtained have negative real parts (i.e., they decay in 
time). For example, the fundamental mode for the Fourier component m is 
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16 2 8 
81° "3 



1/2 

(40) 



This is correct for any a< 4, other values of this parameter being un-physical for the problem at hand. So actually 
the fundamental mode solution is 



v[ (r, z,(/),t)=^r Vi( m ) (r, C)e 



imcp 



S +e s U T+lm * + S_e s (* 



(41) 



where S+ and S- just are (integration) constants. Details can be found in Appendix [Cl 

Clearly, the first order solutions found here lead to solutions exponentially decaying in time . Since in this work we 
arc interested to find transiently growing solutions, we can reasonably make the choice of initial conditions such that 
v[ = p[ = W[ = P[ — at all times (as it was assumed a priori in UNRS). In any case, the decaying solutions, even if 
not chosen to be zero initially, become quickly negligible when the algebraically transiently growing solution appear 
(see next order for a discussion). Moreover we stress that, whatever choices are made, they do not affect at all the 
solution at the next order fsee [jSl [471 H51 and H9")) . 

D. Order e 2 

The second order equations are complicated (see Appendix IB 3|) , but taking only the " unperturbed" parts (base 
flow), which are axisymmetric and steady, one gets the following equation set: 



= - 



- 2rO,^9 = - 



r 2 dr 
8W 2 3z 3 
dz 2r 5 



dW 3z 2 
~dr~ + 2^ 4 

— — (r 3 

r 3 po dr I ^ dr 



Po dz 

dn 



dui 
dz 



1 d 



Po dz \ dz 



dn. 



4 d ( dv2_\ 2_d_ 

3po dz \ dz ) 3rpo dz 



d(ru\ 
] ~dr~ 



1 d 



rpo dr \ ^ dz 



dui 



(42) 
(43) 
(44) 
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= -7T- (tpqu-l) 
r or 



d(ppV2) 
dz 



(45) 



The above equations (|42j|45p are identical to the corresponding steady equation set of KK and UNRS (22-25). 

After subtracting from the complete set we are left with the equations for the non-axisymmctric time-dependent 
perturbations 



~dT ' 


h d<t> 


= 2rtt Q 2 


on' 2 




ui d 
r 2 dr 


dv' 2 
dt 


dv ' 2 


dW 2 

dz 


dp'2 
dt 


dp ' 2 


--!-< 

r or 



J_d_ 

Po dz 

4 d 
3/9o dz 

2 d 
3rpo dz 



du\ 
dz 

i_d_ 

Po dz 
dv'. 



dz 



dtt' 2 
3po <9z 



i__d_ 



9(ru^) 



1 d 



rpo dr V ' <9z 



5^; 



Note that the LHS of all the four equations again contains the operator V^t = dt + flod^,. 



(46) 
(47) 

(48) 
(49) 



Reformulation in terms of linear differential operators 

(|49[) . as in the axisymmetric case, 



We notice that Equations (j4"6"l) - (|4"7]) dynamically decouple from Equations 
treated in UNRS. The whole set of equations can thus be be reduced to 

Vu[ = 
VQ 2 = 



(50) 
(51) 
(52) 



The linear differential operators V.L^J-^ are identical to those in UNRS, in particular the operator C is the one 
inferred (in the previous order) from equation (|35p but here we have, in addition, the non-axisymmetric operators J 
and 7i. The operators are given, in terms of the known function po{z), r\ (a known function of po, r and constants) 
as follows 



V 

C 
T 
Q 

n 

j 



l d f d 

V^t -z- Vtt- 

po dz \ dz 



u 4>t 



4 1 



-V 



d ( d 



3 po t dz 

2 I d d 

3 rpo dz ^ dr 



2 d W d 
oo 

2 1 d 

3 po dz 
5d_ 

3dz 



3 dz rpo dr 

d 

v-. 



2/3 d 
Po *1 



dz 
1 d 



3 
d 

rpo dr oz 



5 2/3 d 

-Po 



2 

dz 2 



25 d( P y z ) d 

6 dz dz 



hd\^_ 

2 dz 2 



>l 



d 2 



d(j> J po d<j>dz 



(53) 

(54) 
(55) 
(56) 
(57) 



where we have used the expressions for the n = 3/2 polytrope, and r\ is as given in eq. ([38 



2. Solution 



We first find the eigenfunctions of the operator V (as in UNRS 50-51). We then decide to consider just the 
fundamental modes to determine u[ and tt' 2 . These can be substituted in the inhomogeneous equation (f5"2"]) to get 
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v' 2 , in view of the the similarity between J 7 , Q in UNRS and 7i,J here (the analogous solution was called in UNRS 
driven acoustics) . We shall discuss here only the case of (the details can be found in Appendix [D]) , since the case 
of n' 2 is the same . 

Mindful of the 2ir— periodicity in <f> we use the Ansatz 



u[(r,z,(j),t) = ^2 ui {m) {r,z)e p 



T+im4 



C.C. 



(58) 



and 



tt 2 (r,z,<f>,t) = ^2 &2(,n)(r,z)e p 



T+im<p 



C.C. 



(59) 



Note that the eigenvalue p is different from s, while T = i/r 3 / 2 is defined as before. 

The substitution of the Ansatz into Vu\ = gives rise, for each m, to a differential equation that is somewhat more 
complicated than the Gegenbauer equation we had in the first order: 



1 f) 2 



, d 3, 
a ^dC ~ ~2~ {j>+im > 



2 



+ ^ } u 1{m) (r,C) = 0. 



(60) 



Still the equation can be tackled analytically using the series expansion (Frobenius) method. As shown in Appendix 
[D]the spatial part of the fundamental mode (k = 0) has a rather simple structure 



and, similarly 

The appropriate eigenvalue p, following from the surface boundary conditions, can take two values 



So equation ([58]) should be replaced by 

oo 



P+e p U T 



P_ e *V) T 



(61) 



(62) 



(63) 



(64) 



where P± are constants. The same holds for il' 2 We note that the homogeneous part of the linear inhomogencous 
equation f|52[) is identical to the equation considered in the first order (i.e., involving the operator C) and that the 
RHS (the inhomogeneous part) can be found using the solutions for u' x and £1' 2 that we just discussed. As is well 
known, the general solution of the inhomogeneous equation (v' 2 ) is the sum of the general solution to the homogeneous 
equation (denoted by v' h ) and a particular solution to the inhomogeneous one (v' p ) 



v 2 = v 'h + V ' P - (65) 

We already know that v' h is an exponentially time decaying function. Since v 2 has to be 27r— periodic in <p we make 
a Fourier expansion. Due to the linearity of the operators, v'^v' ,u'i and Cl' 2 must have the same 4>— dependence. The 
structure of the operators C,J-,G,Ti and J points to a particular solution in the form 



oo oo 
k=0 m= — oo 



\^(r,Qe p ^ +V$ ) + {rX)Te p ^ 1 +v^{r,Q^ + V^~(r, ()Te p (™) x I +c.c, (66) 



(fc)- 



with T and ( as before and where we include, for generality, also the overtones k 7^ 0. Inserting the above expression 
for v' p in equation (|5"2")) and substituting on the RHS the solutions based on equations l|58p. (fBT))) . (|6"Tj) and (|6"2"|) , we 
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are led to a rather long and involved expression (see Appendix [Dl equation ID6[) . The fundamental mode (k = and 
dropping the vertical mode superscript) has the following form (the superscript ± corresponds to the two eigenvalues 



(m) 3(m) 



MC 3 



l l(m) 



(r)C, 



(m) 3(m) ^ J ^ l(m)^ 



By setting to zero the coefficients of t and each power of £ wc get a^ m Jr), a^ m ^(r), b^^Jr) and 6^ m ^(r). 
One may follow a similar procedure for the density perturbation, using 

P2 = Ph + Pp- 

Following the argument outlined in Appendix D. we get that p' h = 0. The Ansatz 



(67) 



(68) 



oo oo 

\ e im4 



fc=0 m= — oo 



c.c. 



(69) 



can be substituted in equation (|49|) to find for the fundamental mode 



Pfm) 



R 



(m) 



= (i-O- 
= (i-O* 



C 4(m)C + C 2(m) 



(r)C 2 + c, 



0(m) 



(>•) 



(OC 4 + 



(70) 



The functions Oj( m ) >^( m ) j c i( m ) an d ^( m ) are known analytical, generally complex, functions of r and depend also 
on the parameter a. Being very complicated expressions, we handle them using the symbolic tools in Wolfram's 
Mathematica 6, but we do not write them in the paper explicitly. Finally we can rewrite, in a more compact way, the 
fundamental driven (by the horizontal velocity perturbations) vertical acoustic modes as 



(m) 



V, 



(m) 



>) 
l(m) 



"3(771) 

bt , 

3(m) 



c 3 , 



(71) 



and 



P%) 

(m) 



= i-C 



2\ 2 



°0(77l) 
U 0(77l) 



U 2(m) 

a 2(77l) 



c 2 + 



C 4(7n) 
U 4(m) 



(72) 



IV. TEMPORAL EVOLUTION OF THE PERTURBATION ENERGY 



Armed with analytical expressions for the dynamical variables of the disk, that are solutions of an IVP in which the 
initial conditions are all small perturbations on a steady base flow, wc may now consider possible physical implications 
of our analysis. 

The solutions we found point to an exponential time decay of most variables, but there is also a transient fast 
algebraic growth (before an ultimate, rather slow, exponential decay) of disturbances of order e 2 in two variables — the 
vertical velocity v' p and density p' p . These quantities can be naturally interpreted as sound waves, which were called 
in UNRS driven acoustics. UNRS study was axisymmetric, while here wc have allowed for azimuthal dependence. 
Exploiting the obvious 2it periodicity in </>, however, we took recourse to Fourier series expansions in this angle variable 
and thus were able to consider separately the various Fourier modes. 

As in UNRS, we perceive the acoustic energy as a relevant variable to follow, but here we consider the time 
dependence of the energy contained in various Fourier modes as well as the effect of different values of a and of 
the (free) radial form of the initial conditions. Following Rayleigh's book [3^| (Chapter XI) we define the acoustic 
energy density as the sum of the kinetic and potential energy densities. The potential term is the work gained/lost 
during expansion/compression and can be easily found using P'dV and P' = p^c 2 sn p' . Hence, the acoustic energy 
volume-density in the m-th Fourier component is 



1 



Pov 2 {m) 



1 r 2 7? 

1 C sOP(rn) 

2 po 



(73) 
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The total acoustic energy volume-density is 



e(r,z,(p,t,(j);a) = -p v + 



2 Po 



(74) 



and we shall consider the fundamental k = mode only. 

The zeroth order vertical velocity is zero, as well as Vi(r, z, <f>) and pi(r, z, </>). Thus remembering that / = f(r, z, 4>) + 
f'(r, z, cj), t), v' 2 = v' h + v' p and p' 2 = p' p , we get 

i = ev[+e 2 (v 2 + v' h +v' p ), p = ep[ + e 2 (p 2 + p' p ), (75) 
where (we repeat all solutions here, for completeness) 



C.C. 



m— — oo 



P'l 



E e !Tn0 ^ (m) (r, C) [^e<-' T + S h _e s ^ T 

,= — ( 
CX) 

E 

,=— ( 

oo 

E 



m — — oo 



m— — oo 



E 



im4> 



*Uh ^ ePim)T + V U r ^ OTe^ T+ + v- m) (r, C)eV.» T + V^r, QTe^r 
^ 1(m) (r,C) [z + e s U T + Z_e s ^ T ] + c.c. 



+ c.c. 



where S±,S^,Z± are (integration) constants, and where 
T = n t = t/r 3/2 , sf m) = --a-im±i 



16 2 8 

— a 

81 3 



1/2 



P, > = Q — ITO ± 1. 



(76) 



(77) 



We recall that 



\h)m are the appropriate Gegenbauer polynomials, while Vr m y V{ m ) are oc ^ polynomials of 
£ and p^ m y R^rn) an< ^ Pi{m) are even polynomials of (. The terms and v' h exponentially decay while v' p ~ Te pT 
grows linearly, so that the algebraically growing terms (like ~ T) dominate the particular solution. Moreover, when 
its eventual decay takes over, the former decaying terms are negligible (for the relevant very small values of a). So, 
for simplicity we have decided to take the initial conditions such that v[ and v' h are zero and thus remain so for all 
times. The same reasoning is valid for p. Therefore when we calculate the acoustic energy we can take just v = v' p 
and p = p p . 

In order to study the temporal behavior of the acoustic energy, we define two integral quantities the radius- 
dependent — i.e, averaged over a ring — energy per unit area of the disk £ r and the total acoustic energy of the 
fluctuation E„ 



fh(r) />27r 

£ r (r,t; a, to) = / / £(r,z,t,<j);a,m)dzd4>, E a (t;a,m)= j £ r r dr. 

-h(r) JO 



(78) 



A. The azimuthally averaged energy surface density 

We obtain £ r for a particular Fourier mode to in the form 



£ r (r,T;a,m) = e~- al F(r\ a, to, cos 2T, sin 2T), 



(79) 



where the function F is known analytically, but the expression is extremely long and we shall not write it out here 
explicitly The derivation is outlined in Appendix F. Instead, we display the results graphically. 
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1 10 100 1000 

T=t a 



FIG. 1: Time evolution of the surface-density of the acoustic energy in the fundamental k = mode at r = 1, a = 0.001 and 
A(r) = C(r) = e' w/4 . The three analytical curves are for m — (solid lower black line), m — 1 (black dotted line), m — 10 
(solid upper red line) and are shown in a log-log plot. £ r {T) is scaled to its corresponding value at T = 0. Notice that £ r is 
modulated by fast oscillations. 




FIG. 2: Same as in Figure [T] but for a single Fourier component (m = 3) and different values of a . Four curves, for 
a = 0.0001, 0.001, 0.01, 0.1 (from the top to the bottom), are shown. 

In the first three figures the temporal behavior of the normalized acoustic energy surface-density contained in 
a particular Fourier component, £ r (r, t;m, a)/£ r (r, 0;m, a), is shown. This quantity is shown as a function of the 
similarity variable T, defined before and for the fundamental k = mode, in a log-log plot. In all the figures the free 
radial functions are set to A(r) = C(r) = e I7r / 4 , which means that the initial perturbation is taken, for simplicity, to 
be r-independent. 

In Figure [TJ a = 0.001 and r = 1 are fixed and the various curves are for different Fourier components. The most 
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t 



FIG. 3: Same as in Figure [2j but for one fixed a — 0.001 and at different radial positions. Four curves for r — 1, 3, 5, 10 (from 
the top to the bottom) are shown. 

significant feature of the behavior displayed in this figure (as well as in all subsequent ones) is the strong transient 
growth (by orders of magnitude) of the perturbation energy, which typically takes place for a rather significant time. 
The non-axisymmctric modes display maximal growth. All the mode have fast oscillations, atop the much slower, 
secular growth. From this Figure it appears that the maximum is attained at t ks 1000 (since r = 1, we have here 
t = T). 

Figure [2] shows the dependence of the transient growth on the value of the a parameter, for the single Fourier 
component m = 3 (other modes display a similar behavior). The radius is fixed, as before, at r = 1. We clearly see 
that the lower is the value of a, the higher is the maximum, and the later it occurs. While a = 0.1 does not give 
rise to any growth at all, for a = 10~ 4 the growth is enormous — by a factor of a few times 10 5 . Since the growth 
occurs when taking into account the e 2 terms of our expansion and the energy is composed of squares of these terms, 
the validity of the expansion is marginal for such a growth (assuming e ~ 10 -2 ). However for a ~ 10~ 3 (close to a 
"realistic" value, as found numerically in the sub-critical hydrodynamic transition (l3j , the growth is somewhat less 
than 10 4 and the asymptotic expansion reasonably holds. 

In Figure [3] a = 0.001 and the Fourier component m = 3 are fixed, while we show the behavior at different radii 
(from the top to the bottom r = 1,3, 5, 10). The growth is higher and it reaches its maximum earlier at smaller radii: 
clearly at large radii the perturbation becomes negligible. 

Before moving on to the description of the behavior of the total acoustic energy, we would like to remark that 
although our choice of the radial functions A(r) and C(r) may seem non-physical, because it is difficult to imagine 
an r-independent perturbation, it is still meaningful. Indeed, since so far we have dealt with the angle averaged 
r-dependent surface density of the energy, it is obvious what to expect when an r-localized perturbation is considered 
instead. The results can be simply scaled, depending on the relative value of the perturbation at the particular radius 
where they are sought (e.g., see Figure[3]). In the next subsection, when we shall consider the total (also r-integrated) 
energy in an extended ring, we will use perturbations that are r-localized (in the form of a Gaussian). 

B. The total acoustic energy in a ring 

The total acoustic energy of the perturbation is evaluated as an integral over r of the afore considered function 

E a (t;a,m)= / £ r (r,t; a,m) r dr. (80) 

-^f mill 

The ring we consider is between r m j n = 1 and r max = 10, far from the inner and outer edge of the disk. 
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3500 




T T 

FIG. 4: Log-log plot of the total acoustic energy in a ring. The left panel is for different modes (m = 1, 3, 10 from the bottom 
to the top). The right panel is for m = 3 and different viscosity (a = 0.0001,0.001,0.001 from the bottom to the top). For 
details see text. 

For simplicity, the radial functions of the perturbation will be taken in the form of Gaussians equal to each other 

A(r) = C(r) = e ™/* e -{r-r f/*-\ ( 81 ) 
so that ro is the center of the perturbation and A its width. 



In principle, our expansion is valid up to a time of the order of e~ 2 . However in making the choice of ro and 
A we should take into account the fact that, in general, the perturbation may propagate with the speed of sound. 
While the boundary conditions on z have been chosen in a physically sound way, we do not have specified any r- 
boundary conditions (this being also the property of the KK and UNRS solutions). Thus, if our ring is determined, we 
should not allow any wave to reach these non-physically natural boundaries, so as not to create spurious effects (e.g., 
reflections that may artificially reinforce the perturbations). Therefore the validity of our results should be limited 
in time. We shall now estimate the time r during which our result is valid. Let Sr be the distance from ro to the 
nearest ring edge: then in our units an estimate of the time for a sound wave to reach that edge is given by 

= e- 1 Sr, (82) 



c s Q (r n )h(r ) 

where eo = M r o)/?~o- Our steady-solution, far enough from the zero torque radius, has h oc r (see the previous section 
and KK) and thus we may substitute eo = e. 

For example, if we take ro = 2 and the inner edge is at r = 1, we have Sr=l and sot~ V2/e . Hence, for e = 10~ 2 , 
the result should be reliable for only a little over 100 time units. The width of the perturbation further limits this. 
A better choice would be to place the perturbation at the center of the ring, i.e., ro = 5.5: then we would have 
t ~ 4.5\/5.5/e, which would give a validity of ~ 1000 time units. 

In Figure |4] we take ro = 4, r max = 10 , A = v2 and plot the evolution of the normalized energy in a ring, 
E a (T)/E a (0). The left panel is for different modes and it shows that the relative growth is higher for higher m. As 
expected, the transient growth is more pronounced for smaller viscosity (right panel). 



V. SURFACE DENSITY SPATIO-TEMPORAL BEHAVIOR 



The formalism developed in UNRS and this paper can be used to follow the time evolution (and, in particular, 
the transient growth) of various small perturbations, which arc included in the initial conditions of the appropriate 
IVP. Wc have already shown the copious transient growth of the acoustic surface energy (averaged over the azimuthal 
angle) as well as of the total acoustic energy in a finite ring of the disk. As suggestive as these results may be, for the 
possible disk energetics, they do not contain explicit spaiio-temporal dynamical information. The possibilities to gain 
the latter are rather abundant, and it would be outside the scope of a single paper to examine a great many of them, 
in detail. Thus, we have decided to conclude here, by calculating and presenting just one of the important dynamical 
variables — the surface density. Other information that can be extracted from our three-dimensional analytical solution, 
found in this paper, will be presented in later works. 
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Expression for the leading terms 



In Section IV we have already estimated that the algebraically growing terms will be the dominant ones. We shall 
repeat now this argument applied to the density. Focusing on the fundamental vertical mode k = (thus dropping 
the superscripts) and taking only one, the m-th, say, Fourier component, we get for the particular solution for the 
density perturbation, appearing in second order in e, the following expression, resulting from equation (|77p . 



P P 



~+ ( \ pt„,\T+ irn <i> 



+ c.c, (83) 



where T = tr 3 ^ 2 , £ = z/h(r) and 



(m) 



-8a/ 5 — i(m =F 1) 



The full asymptotic series for the density includes, however, more terms and is rewritten here, up to the second 
order in e 



p(r, C, <f>, t) = Po (r, C) + e [pi(r, C, <t>) + p'i(r, (, ^ t)] + e 2 [p 2 {r 7 £ <t>) + p' h {r, (, </>, t) + p' p (r, (, </>, t)] 



(84) 



It was shown before that the first order steady term, that is p±, can actually be set to zero. We can approximately 
ignore the steady second order term p 2 j & s well, because even if it is not zero, its magnitude is of the order e 2 as 
compared to the zcroth order steady term. Substituting also the relevant terms in (75) from the explicit formulae for 
Pi and p' p given by formulae (|77p . we get a complicated expression for the density, including the steady base state and 
the mth Fourier mode perturbation. Rather than presenting these expressions, we recall (see the discussion following 
equations ([77]) ). that after some time only the transiently growing terms (those proportional to T) will be dominant. 
These terms will be large for a rather long time, of the order of I /a. Thus, in this example we shall only examine the 
spatio-temporal behavior of these terms. Thus we shall consider 



p(r, C, 4>, t) = p (r, C) + e 2 [Te*™* [A+ m) (r, Q)e^U T + R~ m) (r, Qe^ 



c.c 



}• 



that is, 



p(r,(,(t>,t) = po{r,Q 



Texp ( -g«T 



K(m)( r ' C)e- l[(m - 1)T -"^ + £ ( - m) (r, Q e -i[(^)T-m 4 



c.c. 



To obtain the surface density one has to integrate over the disk thickness, thus 
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E(r,0,t) = E (r) 



Texp 



-aT 



(85) 



(86) 



(87) 



with TZf n {r) = J\ Rf m) (r, ()d( and 6f n = (m T 1)T - m</>. 

Using the expression for po, as in equation (14) of the paper, we can easily get 

E„(r) = h- 3 ' 2 h\- 9 ' 2 f 1 (1 - C 2 ) 3/2 d( w 0.1h 3 (ry 



9/2 



where h ks h\r (see UNRS) and h\ is a constant of order unity, very weakly dependent on the mass transfer rate and 
a. This is a particularly good approximation for r 3> r+ (as we assume). In the example considered here, we shall 
take hi = 1 for the sake of simplicity (the mass transfer rate can always be chosen accordingly). Thus we have 



SoM ~0.1r- 3 / 2 , 

so that the surface density of the unperturbed disk increases significantly for small values of r. 
Using equation (|T2j) we get 



(89) 



(90) 



where, as explained before, the radial function are very complicated, albeit analytically known expressions. Integrating 
over £ we get 



K±(r,a)= / Rf m) (r,0d( = ir 



1 



: 0(m)( r ) + 7 d 2(m)( r ) + o d 4(rn)( r ) 



(91) 
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where we have explicitly reminded the a dependence of this quantity. 

Thus the total surface density (base flow + perturbation) after a sufficiently long time, when the algebraic term 
dominates (i.e., for tj> 100 time units) and before the overall exponential decay takes over (i.e., for t^ 1/a) can be 
well approximated by 

£(r,0,i) = 0.1r- 3 / 2 + 2e 2 |tcx P [~oir \ [cos0-3t(ft+ ) + sin0 m S(ft+) + cos 9+MK m ) + sm0+3(ft-)]} . (92) 

With the definitions of IZ^ and 8^ as above, we can now calculate the spatio-temporal evolution of the surface density 
in the time interval where our approximations hold. 



Results — example of a pattern evolution 

In what follows we shall present graphically three time snapshots of the ratio between the perturbation of the surface 
density (i.e., the e 2 term of the above equation) and the unperturbed surface density, E . It should be remarked that 
the radial functions include, in principle, two functions, A(r) and C(r), which are technically arbitrary, and can be 
only determined by the initial conditions. In the case of the density perturbation only A[r) is needed. We take it to 
be real, for simplicity, and to consist of a Gaussian peak, centered around some radius in the disk, ro, in the region 
we wish to consider. 

A(r)=g er^/ A \ (93) 

where the width A is chosen appropriately. This is done in an effort to mimic a perturbation localized in r. The 
parameter go determines how large the perturbation is at its peak value. Because this expression is included in the 
e 2 term, the resulting initial perturbation is very small, if go is kept to be < 1. 

In this example, we follow a disk ring and are thus far enough away from the zero torque radius, which our scaling 
ensured is much smaller than 1. As mentioned before, we consider an appropriate ring, because both the disk 
unperturbed surface density £o( r ) and the radial function lZ(r) exhibit a power-law decay — the former oc r~ 1,5 and 
the latter oc r~ 2 . Thus the perturbation becomes relatively less important, as compared to the steady quantity, for 
very large radii, while for small radii it may be too big. For the sake of a clear demonstration of the transient growth 
we choose, as before, a = 10~ 3 . 

In the example, for which the surface density evolution is displayed in the figures, we chose only one Fourier mode 
(m = 2) for simplicity. The perturbation was introduced with a real value of A(r), as given above, with go = 0.5 
(other choices of this parameter do not change the result significantly, as long as the parameters remain of order 1) 
and with ro = 5, A = 1. 

As it can be seen in the in Figure fVl which displays the relative surface density perturbation, that is, 

^^^MilMl, (94) 

at time t = 10, in units of fip -1 ^*), the perturbation is still very small. Even though we cannot be sure that at such 
an early time our neglecting of other terms is justified, we are confident that the transiently growing term is still 
insignificant. We also verify that it has the typical m = 2 Fourier component form. 

As mentioned before, we can follow the time-evolution of the surface density perturbation by calculating, with 
the help of Mathcmatica 6, the appropriate analytical expressions. The snapshot figures are also produced by that 
software. The results for t = 100 and t = 200 arc displayed in Figure Ivl and Ivl Two basic features are immediately 
apparent 

1. The absolute value of the perturbation grows with time, as it should according to the non-modal transient 
growth process (see above, in the body of the paper). The growth, at these times, before the exponential decay 
takes over, is approximately algebraic with T. 

2. Since the time variable always appears in the similarity variable T = f2o(r)t = r~ 3 ^ 2 t, the initial perturbation 
pattern (having an m = 2 Fourier angle dependence) is being wound-up by the close to Kcplcrian flow, producing 
successive peaks and troughs in the surface density. 

The winding of the perturbation pattern becomes rather strong for later times, giving rise to a very narrow "radial 
wavelength" of the basic pattern. The overall pattern becomes close to an axially symmetric one by t — 100. Note 
that a relates to the total time-dependent disk surface density according to, 

Z(r,<t>,t) = Z (r)[l + o-{r,ct>,t)}. (95) 
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FIG. 5: The relative perturbation in surface density, a(r,cj>,t), as a function of position, calculated in a ring 1 < r < 7 of the 
disk, at a short time after the initial condition (t = 10). For details see the text. 



(=100 




FIG. 6: The relative perturbation in surface density, a(r,<j),t), as a function of position, calculated in a ring 1 < r < 7 of the 
disk, at time t = 100. The magnitude of the relative perturbation grows and its pattern becomes more complex. For details 
see the text. 



Thus, as long as \a\ < 1, negative densities arc not encountered. Since the pattern at t = 200 is almost axially 
symmetric we can display the total surface density variations along a radial cut, that is to say, on a two dimensional 
plot with respect to radius for a fixed value of the azimuthal angle. In Figure [8] we present such a cut through the disk 
at (f> = 0. The total surface density displays a pattern resembling cylindrically symmetric "waves" , with an amplitude 
of the order of the unperturbed surface density itself (also shown). By t = 200, the variations in the surface density 
are of the same order as the unperturbed surface density, approximately 0.5 EoM^S ^( r > 1.5£oM- 
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FIG. 7: Same as Figure 6 but for time t = 200. 
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FIG. 8: The perturbed surface density (in arbitrary units) as a function of r (solid line). The unperturbed surface density 
So(r) is shown by the dotted line 

VI. SUMMARY AND DISCUSSION— A NEW DIRECTION TOWARD DISK TURBULENCE? 

In this work we considered the approximate nonlinear dynamics of a disturbed hydrodynamical viscous thin disk. 
The base flow is a Keplerian polytropic accretion disk with vertical structure (the KK analytical solution). By means 
of an asymptotic expansion in the small parameter e (the ratio of the characteristic height to radius of the disk) we 
find the temporal evolution of global non-axisymmetric perturbations. While in the first order all the variables decay, 
in the second order the perturbed density and vertical velocity display a strong transient growth. In a short time 
successive peaks and troughs appear in the surface density, similar to what was observed in the axisymmctric study of 
UNRS. The fact that these structures appear for general non-axisymmetric disturbances promotes the conjecture that 
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this phenomenon could be commonplace in non-magnetized ADs and this, in turn, has very interesting consequences 
for dynamics on the small scales. 

To be more concrete, we start by observing that according to (|87|) the surface density has the functional form 
~ Te~ aT sin 6^, where the similarity variable T = tr~ 3 / 2 and the radial "wavelength" 8^ ~ (m =p 1)T. We note 
here two points. Firstly, the amplitude of surface density patterns increase steadily during the algebraic phase of the 
structure's response before viscosity finally gains importance and the disturbance subsequently dies away. The steady 
shortening of the acoustic pattern's radial wavelength is a direct consequence of the non-normal nature of the operator 
governing the pattern's response and the non-separable nature of these solutions. Secondly, the temporal decrease of 
the radial wavelength is enhanced for larger values of m, and we interpret this as the pattern getting "wound up" . 
Thus, we remark that as the pattern winds and the crenelation deepens, the amplitude of the surface density pattern 
grows with time. This latter effect is causally tied to the growing response in the vertical velocity which comes from 
shear energy being converted into vertical mechanical motions. 

The pathways by which this energy is fed into vertical motions can be ascertained by studying the driving terms in 
the inhomogeneous operator (|52|) . Referring to the definitions given in (|54M57|) . we see that the effects responsible for 
bringing about the algebraic time dependence derive from the operators T and Q which represent, respectively, the 
work done by radial compression of the radial velocity u\ and the radial gradient of rz viscous stresses. Technically 
speaking, the operators T and Q 1 each having gradients with respect to r, bring about factors of T(t) when they act 
on the solution to u\ We note, however, that as the viscosity parameter a decreases, the algebraic behavior is 

more strongly tied to the compressional work rather than the work extracted from the viscous stresses, as the latter is 
scaled by a. Nonetheless, all perturbations that include disturbances in the angular velocity and/or radial velocities 
(i.e., u'j^fij) will give rise to this algebraic behavior. However, and although the point may be academic, we remark 
also that the azimuthal compressional work of the angular velocity does not contribute to the algebraic growth. 
This growth specifically comes from the radial compression effect d r u'^. 

The development of structure on finer scales as time goes forth also means that terms in the asymptotic expansions 
formally start to break order. This comes about because the very same radial derivatives discussed above no longer 
remain an order e smaller than the z-derivatives when the algebraic part of the growth becomes substantial and 
this is physically related to the steady development of radial structure [401 ] . This situation becomes especially severe 
for smaller a, as the algebraic growth persists for timescales which are inversely proportional to the viscosity. The 
breakdown of the asymptotic expansion means that the equations of motion must be re-expanded in order to handle 
the evolution of these highly wound structures with fine radial structure. The resultin g eq uations will be something 
like the shearing box equations [U| or some other appropriate model set fill l42l l43lk4 l45j . 

Fortunately, there are preliminary indications of what may happen under these circumstances. There is recent 
literature devoted to exploring what occurs on small disk scales when there is a sizable deviation from a Keplerian 
flow which, in turn, is related to significant radial variations of the surface density. One of the effects we expect to 
happen from the solution scheme employed in this work is that the fourth order correction to the angular velocity, O4, 
will be algebraically forced in T by P2 via the radial pressure gradient. This can be found in the e 4 order expansion 
of @. As the pattern continues to wind and the crenelation deepens, the correction SI4 will algebraically grow with 
T, so strongly that one of several things could possibly happen 

1. Aside from breaking its asymptotic ordering, D.^ would also eventually provide corrections so strong that the 
composite angular rotation profile will satisfy the Rayleigh criterion for axisymmetric instability in many sections 
of the disk. In other words, for a given small scale disk section under examination, the composite rotation 
profile could conceivably develop a 1/r 2 profile (i.e., the inviscid "Rayleigh Line") or steeper. This would, 
presumably, result in the termination of the algebraic growth and replace it with radial transport arising from 
the axisymmetric instability. A transition may occur even before the Rayleigh line is crossed, as Lesur and 
Longaretti [l3j showed that a subcritical transition into a turbulent state does exist and that the turbulent 
activity becomes more vigorous as the Rayleigh line is approached. 

2. Li ct al. [46[ consider the fate of linearized infinitesimal disturbances in a vertically integrated disk model in 
which there are strong radial variations of the surface density. They find that Rossby wave instability occurs 
when the surface density "bumps" are at least twenty percent above the mean. Non-linear simulations show 
that these non-axisymmetric instabilities turn into long-lived vortex trains which transport significant amounts 
of angular momentum [47} . Prior to the results obtained in this current work, a criticism of this proposed 
scenario would have been to say that it is unrealistic to suppose that disks are spotted with seemingly arbitrary 
bumps of surface density. However, the calculation we have performed in this work shows that surface density 
variations are not only common, but they can grow to large amplitude with sufficient time. A reference to 
Figure [5] demonstrates how a small perturbation can develop into sizable fluctuations of the surface density — 
easily meeting the rough twenty percent minimum requirement needed to trigger a Rossby wave instability. 
Furthermore, this transition is reported to occur as the local flow profile approaches the Rayleigh line from the 
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Keplerian state (see Figure 9 of Li et al. [461]). Preliminary calculations performed in a quasi-3D annular model 
of a disk[42j indicate that these trends are robust. We shall detail these results in a forthcoming work. 

We feel that the Rossby wave instability, which is a non-axisymmetric shear instability, is likely to be strong and 
pronounced under these circumstances. Of course, this assertion must be verified by further study. Whatever the 
outcome may be for the smallest scales, the influence of the largest scales upon the smallest ones ought not be ignored 
in the study of ADs. The importance of this could be deeper than previously realized and it may be one of the 
reasons why ADs are such perplexing structures: while there is an obvious separation of scales in a disk, the solutions 
developed here indicate that, if one waits long enough, the dynamics originating on large scales invariably generates 
power on the smallest scales too. Usually this downscale cascade of power is rationalized as happening because of 
nonlinear mode-mode interactions. However, in this case it comes about due to the fundamental inseparability of the 
solutions to the lowest order dynamical response. Nonetheless, might it be incomplete to examine small scale disk 
behavior without representing the dynamical influence precipitating from the large scales? 

If secondary instabilities do develop as a result of one or more of these processes, then we venture to say that 
such radial variations of the surface density could be cither maintained or replenished due to the anomalous activity 
generated by them. This is because such activity could conceivably generate power back onto the largest scales and, 
consequently, become the seed disturbances for the large scale dynamics elucidated in this work. This would complete 
a dynamical cycle describing sustained disk activity. Even if such small-to-large-scale causal connections are either 
absent or insignificant in disks, one can also envisage that bursts of activity and transport can occur in some disks 
simply due to random perturbations from outside in the way discussed in Ioannou and Kakouris Disks certainly 
do not sit in isolation and periodic disturbances of them by stars passing nearby them is likely, especially for disks 
found in crowded environments like young star- forming regions. 

The outcome of our work can be summarized as follows: 

• Although hydrodynamical thin accretion disks arc linearly stable, we find that the transient dynamics of initial 
3D non-axisymmctric perturbations can give rise to substantial growth. This confirms that this transient growth 
is not restricted to axisymmetric disturbances like investigated in UNRS. 

• In particular, every perturbation that disturbs the radial velocity, leads to an evolutionary phase in which there 
is algebraic growth of the density and vertical velocity. After a longer time this temporal response gives way 
to an exponentially decaying phase wherein viscosity dominates and the perturbed quantities eventually go to 
zero. 

• Due to the non-normal nature of the linear operators involved, the evolution of the perturbation patterns are 
controlled by a similarity variable T, a non-separable combination of the radial coordinate and time. This 
leads to a winding of the perturbations, producing successive peaks and troughs. This directly contributes to 
the resulting algebraic growth of the other quantities discussed above as comprcssional work converts energy 
contained in the wound pattern into vertical mechanical motions. 

• For a given particular Fourier mode, the acoustic energy associated with such a perturbation grows more for 
higher azimuthal number m and for smaller viscosity a. 

• We conjecture that as the perturbed surface density continues to wind up and deepen, secondary instabilities 
could arise (e.g., Rayleigh or Rossby wave instabilities). The development of such instabilities and their interplay 
with the large scale perturbed disk structures could lead to sustained turbulent activity. These processes could 
contribute to the enhanced transport of angular momentum needed to match observations of ADs. 

The aforementioned results should be ubiquitous for general non-magnetized accretion disks. Therefore they are 
relevant for the understanding of accretion disks around compact objects, whether they be active galactic nuclei or 
close binary systems, as well as for protoplanetary disks and circumstellar disks around Be stars. In the future we 
plan to investigate the possible development of secondary instabilities and the observational consequences of our 
findings in the different contexts. 
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APPENDIX A: THE VANISHING OF fii AND u 



Rewriting the first order equations (|27|) and (|28[) with the definitions U = uq and V = 2rf2i gives 



2a 

T 

2a 

T 



3 2/3 



3 2/3 



Substituting now the zeroth order solution p 2 / 3 = (ft 2 — z 2 )/(5r 3 ) and rearranging leads to 

5z 



2 <A 5C/ „ n 

z— + v = o 

<?z 



d 2 E7 

dz 2 V 3 / dz 



ft 2 



52 



ft 2 



V. 



a 



h 2 -z 2 
a 



ft 2 



V = 
[7 = 0, 



(Al) 
(A2) 

(A3) 
(A4) 



where the subscripts z denote here the differentiation with respect to z and a = 15/(2a) is a constant. 

Using now Q{z) = (h 2 — z 2 ) 5 / 2 as an integrating factor for the first two terms in both of the above equations, we see 

that 



d f„dU 







dz 

d_ f„dV 

dz 



-- -a(h 2 - z 2 f' 2 V 



(A5) 
(A6) 



Multiplying the first equation by U, the second by V, adding and integrating over the domain [— ft, ft], gives, after 
dropping the integrated parts, 



dU 

dz 



dV 

dz 



Q{z)dz = 0. 



(A7) 



Because Q{z) ^ 0, except at z = ±h and the functions U, V are bound, they must be equal to constants. Thus, it 
follows from equations (|A3IIA4[) . that U = V = 0, except perhaps at z = ±ft. However, since they are bound and 
constant in all the domain, they (and hence uq and fix) must be zero identically. 



APPENDIX B: COMPLETE 0(1), 0(e) AND (e 2 ) SYSTEMS 

The complete equation sets in these orders will be given here. Complete means here that it is only assumed that the 
zeroth order functions and the non-perturbed portions of higher-order functions are time independent and nothing is 
assumed about the axisymmetry of the solutions. 

We shall also use here the notation for the full function at z-th order, /; = /j(r, z, <f>) + f[{r, z, cf>, t). 



1. Zeroth order system 

The complete equations at this order are relatively simple because there is no time dependence at this order. 

n 2 = 1 (Bi) 
dn Q on d ( dci \ 



dz d(j) dz \ dz 

dv , _ dv z dW 4 1 d ( dv \ ldf dQ \ 21 d ( dQ \ 

v — + n — = — =■ - — — + - — -5- r)— + — — t]— — - — — t]— — (B3) 



dz d<t> r 3 dz 3 po dz \ dz J po d(f> \ dz J 3 po dz 

d(p v ) d(p n ) 



dz 



= (B4) 
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2. First order system 



8u du ft 1 d ( du 

dz d(p p dz \ dz 



(B5) 
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3. Second order system 
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APPENDIX C: FIRST ORDER SOLUTION, DETAILS 



Solution 



Substituting 



v[ (r, z,^t)=Y J «!(«») (r, C)e sO0t+ ™^ + c.c. 



(CI) 



i.e., equation (|39p in the linear differential equation (|36[) . gives the ordinary differential equation 



(i-C 2 ) 



d 2 



5C 



8 



1 + s 



"l(m)(^C) = 0, 



(C2) 
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where s* = s + im. The previous equation (which is similar to B.3 in UNRS, who derived it for m = 0) is valid for any 
to and n = 3/2. We can find the solution with Wolfram's Mathematica 6 (or see [37], [HI)- Using the initial condition 
vi( m )(r,0) = 0, one gets 

v 1(m) (r, C) = A(r \ /4 ■ k 3/2 (0 QT (0) - QT (0 PT (0)1 , (C3) 

(i - C 2 ) 

where P„{() and are the associated Legendre functions of the first and second kind, respectively. A(r) is the 
constant (in Q of integration and is free, depending on the initial conditions. The parameter v relates to the parameters 
of the system via, 



u(a,s*) = (^-V22.5 + 12as* + V 90 - 270^2 + 192as»J / (2V22.5 + 12as*) . 

In addition, in order for the solutions above not to be singular at the boundaries £ = ±1, v must be equal to half 
integers (see next section). The general solution to (|34p is the sum of the general solution to the homogeneous part 
and a particular solution to the inhomogeneous. The homogeneous solution is simply a wave — i.e., any function in 
the form (((j> — Clot). However, if we consider the complete system (|33fl34|l . then we get that — Clot) must be zero. 
Therefore p[ is only equal to the particular solution. It can be found by substituting v[ resulting from (|C3|) in (|34p 
and using the Ansatz 

oo 

p' 1 (r,z,^t)= Pi( m) (r,C)e sO0t+ "^ + c.c., (C4) 

Then, for a given m, we find 



m— — oo 



( m ) I y\ 

Pi {r,Q = 



10V5r 3 / 2 s,{-l + C 2 f /4 



(i-C 2 )W 



{2v - l)A(r)h(r) x 



CP* /2 (0Ql /2 (0) - P*l\(OQl /2 (0) + ^ /2 (0) (~CQl /2 (0 + Q 3 JA(0 



(C5) 



-,3/2 



The eigenvalues, s, for any to, will follow from the boundary conditions on £ (see below). 



Vertical boundary conditions 

We assume that the lagrangian pressure perturbation vanishes at z = ±/i(r) 

^ = 0=(^+«-V)P = 0, at z = ±h(r). (C6) 
Using the polytropic equation and the continuity equation, (|C6[) can be simplified to 

PV-v = 0, at z = ±h(r) (C7) 
At this order, since Uo = Oi = 0, it reduces to 

dv' 

P o (r,z)—± = 0, at z = ±h(r). (C8) 

We substitute the Ansatz for v{ and the solution IC3I into (|C8|) and consider the limit for ( — > ±1. Therefore the 
boundary conditions (for n > an integer or half odd integer) resulting from regularity conditions at the disk upper 
and lower edges are satisfied for ^ = 2fc + 5/2, where k > is an integer or k = —3/2. Following the discussion in 
UNRS B.l we can consider the fundamental mode (k=0) and so find the two fundamental eigenvalues (s^)), that lead 
to a typical oscillatory decay. The values of a which we consider in this work (a <C 1) do not allow a non-oscillatory 
decay — see (14TJ1) . The actual solution for the fundamental mode (k = 0) becomes 
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«i(r, 0, 0, i) = £ *!(«.)(»•> ()(S + e s ^™ + S^W + ^) + c.c, (C9) 



and 



P ;(r,z,0,i) = ]T^ m) (r,C)(^e<^ (CIO) 



— oo 

Surface stress conditions 

It is easy to see that the vanishing stress conditions are satisfied on the surface 



lim r/n • Vu — > and lim r/rh ■ Vf2 — > 0, (Cll) 

js — >/i z — >h 



in which n = z to lowest order. 



APPENDIX D: SECOND ORDER SOLUTION— DETAILS 

As detailed in the text. ([60|) holds for the m-th Fourier component of the function u[(r, z,(j>), see (|58|). We rewrite 
it here for convenience. 



1 m ,2, d 2 . d 3 



12 9 



+ i ^ 1(m) (r,C) = 0, (Dl) 



where £ = z/h(r) and p* = p + im is an eigenvalue. 

A similar equation is also satisfied by the Fourier components of £l' 2 (r, z,<f)), see (|59|) . To solve the above equation 
we assume a truncated series expansion, which for a general k > reads 

«iW=^) = E^MC 2(fe+1) - 2j , (D2) 

We consider the fundamental mode k = and the relative eigenvalue p, as above, and obtain 

ui=u? ] =A(r)Q 2 +B{r), (D3) 

where we have dropped the (m) subscript for economy of notation. 
Inserting (|D3[) into (|D1|) leads to 



C 2 (225 (1 + j£) + 144a(5?j* + 4a)) A(r) + (225 (l + pi) B(r) - 24a(5p» + 4a)A(r)) = 0. (D4) 
We set to zero the coefficient of ( 2 to find p*, and therefore p, and the coefficient of £° to find B(r) 



Ptn)=~l a ~ im± ^ J$j = ~i (D5) 

These modes are decaying oscillations. We then substitute v', z = (,h,U\ as in (|D3[) and the fundamental mode 
2 = ^2 0) = Q(r){z 2 /h 2 - 1/6) in (£2) we get (n = 3/2) 
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J 49Ah 217(p + im)aAh ( iaAh 7. L ^ 2 ., . . , \ 7, 2 , . . , ,, 

+ C[2A/i' + -(p + im)QL4ft / j (D6) 



Vertical boundary Conditions 

At the second order (|C7|) becomes 

Pn / / du\ . „ ,dv!> dfln. „,dv\ „ , , , , ._ . 

— +Po b^ + ^fR^T—O, at z = ±A(r). D7 

r or oz oq> oz 

(2 2 \ 5/2 
h 5 ~i ) is zero at z = ±/i and W2 and Q' 2 and their 

derivatives are finite. The last term vanishes as well, for k = —3/2 or integer > 0. Therefore the BC is satisfied also 
at this order. 

APPENDIX E: MASS ACCRETION RATE 

We want 

p2tt rh 

r I d(f) pudz = —M = const. (El) 

JO J-h 

If we expand p and u, we get different equations at different orders. The time-independent parts are 

/>27T ph />27T ph 

i d(f> rpouidz = —Mi = const. / d<f> / rp 2 uidz = —M 2 = const. (E2) 

JO J-h Jo J-h 

Therefore M = eM x + e 3 M 2 . 

For the time-dependent part we get 

2tt ph p2ir ph 

d<t> I rpouidz = 0, / dip I rp\u\dz = 0. (E3) 

J-h Jo J-h 



/>2tt />h />2tt />h 

/ d(f> rp 2 U\dz = 0, / d(f> / rp 2 u Y dz = 0. 
Jo J-h Jo J-h 



(E4) 
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However 



rp Y u x dz 7^ 0, / d(j> i rp 2 u x dz ^ 0. (E5) 

J-h JO J-h 

Indeed in both cases the function that has to be integrated is even in z and is not 2~k— periodic in (f>. We can satisfy 
the first equation by setting p[ = 0, but the second equation does not vanish and gives a correction to M that grows 
and subsequently decays like ~ e 3 Te~~ aT . This is a fluctuation of the order of e 3 on a quantity (M) of the order of 
e. It can be neglected for e sufficiently small, that is for a not too small. 

APPENDIX F: ENERGY-DETAILS 

The particular solutions v' p and p' p arc both in the form f = A m e Tm +B m e m +c.c, see (|66|) . Let now A m = A R +iAj 
and, likewise, T m = T R + iTj, B m = B R + iB h G m = G R + iGi, where ^ R ,^i,r R , r I; B R ,B 1 ,G R and Gi are all 
real quantities. Then we have 



where 



Re{f p Y = e 2 '«{T a +T b + T c ), (Fl) 

T a = 2A 2 + 2A 2 R + 2B 2 + 2B 2 + (AA^ + AA R B R ) cos (2T) + {-AA^ + AA R B R ) cos + Gi), 
T h = -1(B\ + B R ) cos (2Gi) + (-2A 2 + 2A 2 R ) cos (2Ti) + i(A R Bi - AA Y B R ) sin (2T), 

T c = -A(A R Bj + AjB n ) sin (Tj + Gj) - iB T B R sin (2Gi) - 4AjA n sin (2^). 

v' p and p'p have the same phases, i.e., T R = G R = — (8/5)aT, T\ = T(l — to) + m0 and T\ = T(l + m) + to0. The 
integral of the quantity in (|F2j) can be simplified using trigonometric relations 

/ -Rc(f p ) 2 d<j> = 2ire-^ aT [A 2 R + Aj + B R + B 2 + 2(A 1 B I + A R B R ) cos (2T) + +2(A R B I + AiB R ) sin (2T)] , 
Jo 1 

(F2) 

For ui — the above equation reduces to 

/ iRc(/;) 2 # = 46-^ 7T [(Ar + S R ) cos(T) - (Ai - Bi) sin(T)] 2 (F3) 
Jo * 

We now write equation ([73]) in the form (Fl) and proceed with the vertical integration, remembering that in 
the fundamental mode (k = 0) we had for the radial and angular perturbations Ui( m )(r, C) = ^( r ) (C 2 — |) an d 
^2(m)( r 7C) = (C 2 — g); with the radial functions A(r) and C(r) free. 

Moreover, we notice that h(r) — > (2A) 1 / 6 r = cir for r >> r*, with ci depending on a and on the mass flux. However 
it multiplies every coefficient and therefore we can set it to one without loss of generality. After the integration in the 
vertical direction, we obtain £ r in the form 

£ r (r,T;a,m) = e~ ^ aT F(r; a, m, cos 2T,sin2T). (F4) 
F is a known analytical function: we shall not write it out explicitly for space considerations. 
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